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Convergence of a Fast Explicit Operator Splitting 
Method for the Molecular Beam Epitaxy Model * 

Xiao Li^ Zhonghua Qiao^ Hui Zhang§ 


Abstract 

A fast explicit operator splitting (FEOS) method for the molecular beam 
epitaxy model has been presented in [Cheng, et ah, Fast and stable explicit 
operator splitting methods for phase-field models, J. Comput. Phys., sub¬ 
mitted]. The original problem is split into linear and nonlinear subproblems. 

For the linear part, the pseudo-spectral method is adopted; for the nonlinear 
part, a 33-point difference scheme is constructed. Here, we give a compact 
center-difference scheme involving fewer points for the nonlinear subprob¬ 
lem. Besides, we analyze the convergence rate of the algorithm. The global 
error order 0(r^-|-h^) in discrete L^-norm is proved theoretically and verified 
numerically. Some numerical experiments show the robustness of the algo¬ 
rithm for small coefficients of the fourth-order term for the one-dimensional 
case. Besides, coarsening dynamics are simulated in large domains and the 
1/3 power laws are observed for the two-dimensional case. 

Key words: molecular beam epitaxy, fast explicit operator splitting, finite dif¬ 
ference method, pseudo-spectral method, stability, convergence. 

1 Introduction 

Recently, the molecular beam epitaxy (MBE) has become an important technique 
for the growth of thin Elms. By using the MBE technique, it is possible to grow 
high-quality crystalline materials and form structures with high precision in the 
vertical direction [9]. There has been a large amount of research interest in the 
dynamics of the MBE growth. Different kinds of models have been developed 
to describe the growth evolution, including atomistic models, continuum models, 
and hybrid models [6]. In our work, we are interested in the continuum models 
for the evolution of the MBE growth. The evolution is governed by the following 
nonlinear partial differential equation: 

ut = V ■ [(|Vm|^ — l)Vn] — (5A^n, {x,y) e D, t G (0,T], (1.1) 
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where 5 > 0 is a constant, hi = (0, 2L)^ with L > 0, and m : hi x (0, cx)) —)■ M is an 
r2-periodic scaled height function equipped with the initial data 

u{x,y,Q) =Uo{x,y), {x,y)eVL. 

The fourth-order term models the surface diffusion and the nonlinear second-order 
term models the Ehrlich-Schwoebel effect [5, 15, 20]. The equation (1.1) is the 
gradient flow with respect to the L‘^{VL) inner product of the energy functional 

E{u)= f QdVwp - 1)^ + dxd?/. (1.2) 

With the periodic boundary condition, it is easy to show that the energy E de¬ 
creases with respect to the time. For the coarsening dynamics governed by (1.1), 
the exponents measured experimentally are 1/3, which is observed in numerical 
simulations of the MBE growth [1, 19, 23, 29] and analyzed quantitively by intro¬ 
ducing a kinetic scaling theory [14]. 

There have been many theoretical and numerical studies on the MBE models. 
The well-posedness and regularity of the initial-boundary-value problem of the 
model (1.1) are studied in [14] using the Galerkin approximation method. For the 
MBE simulations, a large computational domain is necessary in order to make 
the effect of periodicity assumption as small as possible and to collect enough 
statistical information such as mean surface height and width of the pyramid-like 
structures. Besides, a sufficiently long integration time is necessary in order to 
detect the epitaxy growth behaviors and to reach the physical scaling regime. To 
carry out numerical simulations with large time and large computational domain, 
highly stable and accurate numerical schemes are required. The equation (1.1) 
is highly nonlinear with a small surface diffusion parameter S, which makes it 
difficult to design an effective numerical scheme. In [17], two stable and convergent 
linearized difference schemes are derived by using the method of reduction of order 
[25]. The convergence rates are 0{t + and 0{t‘^ -|- h^) in discrete L^-norm, 
respectively. Both the nonlinear part and the diffusion term are treated explicitly 
there. In [19], two unconditionally energy stable difference schemes are presented. 
These two schemes are second-order convergent in time and nonlinear. Because of 
the unconditional stability, an adaptive time-stepping strategy is purposed there. 
In [23, 27], the Erst- and second-order (in time) convex splitting schemes are 
constructed under the framework exploited by Eyre [7]. Still, both the two schemes 
are nonlinear and unconditionally energy stable. The similar technique has been 
used extensively on different phase held models, e.g., the phase field crystal model 
[28], a diffusive interface model with Peng-Robinson equation of state [16], etc. 
In [29], the authors introduce an implicit-explicit scheme combined with Fourier 
pseudo-spectral approach, where the nonlinear term is treated explicitly and the 
fourth-order term implicitly. To guarantee the stability, they add an extra artificial 
term consistent with the truncated errors in time. However, the condition, under 
which the energy stability can be obtained without any restriction on time step, 
depends on the unknown numerical solutions. In [18], a mixed finite element 
method with Crank-Nicolson time-stepping scheme is presented and the energy 
laws are proved for both semi- and fully-discrete form of the scheme. 

In [1], a fast explicit operator splitting (FEOS) method based on the Strang 
splitting schemes [24] is constructed to simulate the MBE equations for both one- 
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and two-dimensional cases. The main idea of the method is to split the original 
equation (1.1) into nonlinear and linear parts whose exact solution operators are 
denoted by S^f and Sc, and then to evolve one splitting step (from t to t + t) via 
three substeps: 


u{x,y,t + T) = Sc(^^'^S^{T)Sc(^^'^u{x,y,t). 

A similar strategy has also been used to solve the phase held crystal equation in 
[13]. In [1], the nonlinear part is solved by the 33-point center-difference scheme 
combined with the large stability domain explicit Runge-Kutta solver, and the 
linear one is solved by the pseudo-spectral method. Their numerical experiments 
indicate that the proper constant time step should be r = 5/100. In addition, 
the FEOS method has also been successfully utilized on the convection-diffusion 
equations equations [2, 3, 4] and the modihed Buckley-Leverett equations [10]. It 
is capable to achieve a reliable numerical solutions in an efficient manner, that is, 
only few splitting steps are preformed [4]. 

In our work, we concentrate mainly on the convergence analysis of the FEOS 
method for the MBE equation in the two-dimensional case. The main issue, which 
is different from that in [1], consists of three aspects. First, we discretize the 
nonlinear part by a 25-point center-difference scheme in space and the explicit 
strong stability preserving Runge-Kutta method in time, and combine the so- 
called “frozen coefficient” technique with the Fourier analysis method to derive a 
constraint on the time step for the stability. Second, we analyze the convergence 
of the entire algorithm. The global discrete L^-error consists of the truncation 
errors from the splitting, the nonlinear and linear schemes, respectively. Third, 
we carry out some numerical experiments to verify the convergence rate, and 
test the robustness of the algorithm with small S in the one-dimensional case. 
Numerical experiments suggest that the time step can be set as r = 5/10 using 
our algorithm. This result is a little better than that in [1], because the difference 
scheme for nonlinear part involves fewer points, which may loosen the restriction 
on the time step. Besides, we consider the two-dimensional coarsening dynamics 
to observe the —1/3 power law of the energy and the 1/3 power law of the mean 
height. 

The organization of this paper is as follows. In Section 2, we present the fast 
explicit operator splitting method for the two-dimensional MBE equation, and 
give a sufficient condition for the stability of the algorithm here. In Section 3, the 
discrete L^-error estimate of the FEOS method is shown both theoretically and 
numerically. Further numerical experiments are carried out and the power law for 
the coarsening dynamics is observed in Section 4. Some concluding remarks are 
given in Section 5. 

2 Fast explicit operator splitting method 

Here we present the algorithm developed in [1] where the nonlinear and linear parts 
are approximated by different methods, and construct a more compact difference 
scheme for the nonlinear part. 
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2.1 Splitting strategy 

In [1], the equation (1.1) is split into the nonlinear part 


( 2 . 1 ) 


ut = 'V ■ (|Vm|^Vm), 


and the linear part 


Ut = —Au — 5 A‘^u, 


( 2 . 2 ) 


whose exact solution operators are denoted by Sj^ and Sc, respectively. Introduc¬ 
ing a splitting time step r, the solution of the equation (1.1) is resolved from t to 
t + T via the Strang splitting method [24] consisting of three substeps: 


/T\ 

u(x, 9 ,i + t) = 5cI -j5Ar(’')5 £2 j“(^, V, «)■ 


(2.3) 


In general, if all the solutions involved in the three-step splitting scheme (2.3) are 
smooth, the operator splitting method is second-order accurate [24]. 

For the nonlinear subproblem (2.1), the solution is L^-stable with respect to 
the initial data, which is described precisely by the following proposition. 

Proposition 2.1. Assuming thatuo,vo G we have 

||‘5'A/-(t)Mo - >S'Ar(t)no||L2(o) < ||mo - t'olUaco), Vt > 0, 

where = {u E \ u is -periodic}. 

Proof. Set u{x,y,t) and v{x,y,t) to be the solutions of (2.1) with the initial data 
m(-, •, 0) = uq and v{-, •, 0) = uq, respectively. Let w = u — v, then we have 

wt = V ■ (IVmI^Vm - iVnl^Vn) = ■ ((|Vm| 2 + |Vnp + \Vu + VnnVw). 


Taking the inner-product with w and noting the periodicity, we obtain 


dt 

which leads to 
that is. 


;5l|tt'(«)lll(n) = - / dVtil" + + |Vu + Vt.d)|VTO|'"da:d!, < 0, 


k(f)l|L2(0) < lk(0)||L2(f^), Vf > 0, 


\\u{t) - v{t)\\L 2 {n) < ||m(0) - n(0)||L2(Q), Vf > 0 , 
which completes the proof. 


□ 


In practice, the exact solution operators S^ and Sc are to be replaced by 
their numerical approximations. In the following two subsections, we present the 
numerical methods given in [1], while the algorithm for the nonlinear part is a 
little different. 
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2.2 Center-difference scheme for the equation (2.1) 

Using the method of lines, the nonlinear subproblem (2.1) can be reduced to a 
system of ODEs, which can be efficiently and accurately integrated by a stable 
explicit ODE solver. Here we adopt the fourth-order-difference to discrete the 
space, and choose the third-order strong stability preserving Runge-Kutta (SSP- 
RK3) method [8] as the ODE solver. 

Introducing a spatial scale h = 2L/J, where J = 2N is a positive even integer, 
the grid nodes are dehned as {xj, yk) = (jh, kh), j,k = 1, 2,..., J. The fourth- 
order semi-discrete scheme for (2.1) can be written as [12] 

duj^kjt) _ -Fj+2,j,k{t) + 8T)+ij,fc(t) - 8Fj_ij,fc(t) 
dt I2h 

I —Gj^k,k+2(t) + 8Gj^k,k+l(t) — 8Gj^k,k-l(t) + Gj^k,k-2(t) fo A\ 

' TTTt ? 


where 


Fj-\-iy^k F (JyUx) j^l^j^ki 1 Gj^k,k+i G(y{Ux) j^k,k+ii id^y) j,k,k+(^ i ^ if) 

here F(p, q) = ij? -|- g^)p, G(p, q) = {p^ + q‘^)q, and 

25'Uj-|_2,/c T -I- 3Uj—2^k 


[U 


X )j+2J,k 


{Ux)j+lJ,k — 
id^x) j—l,j,k 

{Ux)j-2,j,k = 


12h 

T Rl^j-t-i,fc T k ^j—2^k 

m ’ 

^j-\-2,k “t“ 18'Ujy lOWj — iy 3Uj—2^k 

m ’ 

~dUjj^2,k T k ‘ddUj^k T dSUj—\ k 25'Uj_2^/c 

V2h ^ 


^y) j,k,k+2 


id^y') j,k,k+l 


\Uy)j^k,k-l — 


^y )j,k,k—2 


I2h 

2dUj^k-\-2 d8xij^k+l T ddUj^k ^dUj^k—1 T dUj^k—2 

m ’ 

dUj^k-\-2 T fO^jy+l T dUj^k—1 ^jyk—2 

V2h ’ 

^y/c-t-2 T 1 dtlj k—2 

V2h ’ 

"3'Uy/c+2 T Ib'Uj'y^l 3dUj^k T dStlj^k—l 2dxij^k—2 

m ■ 


(2.5a) 

(2.5b) 

(2.5c) 

(2.5d) 


f = ±1, ±2, (2.5e) 

+ + ^ ^ ±1_ ±2, (2.5f) 


(2.5g) 

(2.5h) 

(2.5i) 

(2.5i) 


The fully-discrete scheme for (2.1) is obtained by applying the SSP-RK3 method 
[8]. This completes the numerical approximation of the operator Sj^f. We notice 
that our scheme (2.4)-(2.5) is fourth-order in space, which is same as the scheme 
(2.7)-(2.9) given in [1]. In addition, our scheme is more compact than the scheme 
in [1], since the former utilizes 25 points while the latter 33, as proposed in Fig. 
2 . 1 . 

According to the property of strong stability preserving, the stability restriction 
of the SSP-RK3 method is identical to that of the forward Euler scheme. We 
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in [1]. 

Figure 2.1; The “x” represents the point involved in the scheme expanded at the 
“o” point. 


use the “frozen coefficient” strategy to analyze the stability of the forward Euler 
scheme 

_ ,,.n _ T^n _l_ o jpn _ o j^n i jpn 

^j,k _ ^j+2,j,k ' °-^j+l,j,k °-^j-l,j,k + ^j-2,j,k 


12h 

/^n _1_ o/^n _ _L 

~^j,k,k-\-2 “T ^^j,k,k-\-l ^^j,k,k-l ^j,k,k-2 

l2h ’ 

where ^= Fj+ij{tn): Glk,k+e = Gj^k,k+i{tn): i = ±1, ±2, that is, 


+ 


( 2 . 6 ) 




^ = i\^uyy2,,k 

+ iN^\")U,,k- 
+ iWXk,k+2 ■ 
+ iNu\XAk-i ■ 




=^)j+2,j,k 

12h 


+ (|Vn| 


2\n 


2(u 




3h 


[u 


yJj,k,k+2 

V2h 

3h 


+ iNu\X-2,,k 


+ i\^u\Xk,k^i 


+ i\^X)lk,k-2 ■ 


3h 




^F-2,j,k 

12h 

‘^yy)j,k,k+l 


3h 


[u 


yjj,k,k-2 

~12h 


It is observed that the terms {uX+e,j,k, {uy)Xj,k, Mlk,k+e^ (XX+i approx¬ 
imate the values 'Ux{xj-\.^^yk',tn)^ Uy(^Xj-^-£^ykitn}j Uxi^Xj^yk+i^tfi^j Uy(^Xj,yk-\-£jtix) 
with the error OX), respectively. Freezing the prefactors of the square bracket 
terms by the constant 

A = max{(| „ i\Vu\X,k,k+i : i = ±1, ±2}, 

we obtain the following linear scheme: 


y,k 

T 


= A 


-u 


'i+2,fc + ^x+^,k - 3o«i,fc + ^x-i,k - y-2,k 


12/l2 


ulk+2 + 16m? 


yfc+i 


3Xk + 16< 


j,k-l '^j,k-2 


12h^ 


which can be transformed into the following form: 


“ (1 ~ + ~:rX+i,k + X-iM + y,k+i + 


4r 
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“ + '“i,fc+2 + '^],k-2)^ (2-7) 

where r = Ar/K^. Using the Fourier analysis method, the symbol of the difference 
scheme (2.7) is 


T 

p(o'i, 0 - 2 ) = 1 - -[(1 - cos(Ti/i)( 7 - cosaih) + (1 

O 

Therefore, |p((Ti,(T 2 )| < 1 if and only if 


cos a 2 h) (7 — cos < 72 / 1 )]. 


0 < r < 


(1 — COS(7i/i)(7 — COSCTih) + (1 — COS(T2h)(7 — COS(T2h) 


As (1 — c)(7 — c) G [0,16] when c G [—1,1], we obtain r < —, namely, 

16 


r < — -h\ 
- 16A 


( 2 . 8 ) 


Obviously, this is a sufficient and unnecessary condition for the stability of the 
Euler scheme (2.6), and thus, of the scheme (2.4) combined with the SSP-RK3 
solver. 


2.3 Pseudo-spectral method for equation (2.2) 

In [1], the equation (2.2) is solved by the pseudo-spectral method via the following 
procedure. They hrst use the FFT algorithm to compute the discrete Fourier 
coefficients \upq{t)} from the point values {■Uyfc(t)}. Then they calculate Upq{t + 
t) = e^p'>'^Upq{t) , where 

, n\p^ + q^) ,fn\p^ + q^)y 
V- J2 J2 ) ■ 


Finally they recover the point values of the solution at the new time level, {uj^k{t + 
r)}, from the discrete Fourier coefficients {upq(t -l- r)} using the inverse FFT 
algorithm. 

For the self-consistent of our statement, here we give some formulas to be used 
in the next section. For the continuous function u{x,y,t), there exists the Fourier 
series in the complex form at time t\ 


OQ 00 


U[X, 


:.».«)= E E (2.9) 

p=—oo q=—oo 

where the Fourier coefficients are given by 

Upq{t) = ^ [ dxdy, p,q = 0,±1, ±2,.... (2.10) 


'n 


It is easy to see that the Fourier coefficients satisfy the following ODEs: 
^S„(() = A„S„(«), \„= 
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The exact solution is 


Upq{t + t) = e^P'^'^Upqit), p,q = 0, ±1, ±2,..., 


and then 

OO OO 

u{x,y,t + T)= ^ + (2,11) 

p=—oo q=—oo 

It is easy to see that the pseudo-spectral method purposed in [1] is the discrete 
form of the procedure above. 

In the theory of the spectral method [21, 22], the FFT and the inverse FFT 
algorithm can be expressed as 


p,q = -N,...,N, (2.12) 

^ ^p^<i 1=1 k=i 

and 

N N 

“mW = E E (2.13) 

p=—N q=—N 

where Cp and Cq are dehned as 


Cr 


2, |r| = iV, 
1, |r| < N. 


(2.14) 


The psendo-spectral procednre can be expressed as 

u{t + r) = g)}. 


where u(t) is the matrix with the elements {uj^kit) : j,k = 1,2,, J}, and 
are the discrete Fonrier transform and the inverse transform, respectively. 
Using the Parseval’s formula and the fact that < e^s (for any p,q), we 

obtain 

||n(^ +'r)|| < e^||M(t)||, (2.15) 

where || ■ || represents the discrete L^-norm, that is. 


u = 




1=1 fc=l 


The inequality (2.15) implies the stability of the pseudo-spectral procedure. 


3 Error analysis and accuracy tests 

Here we investigate the convergence rate of the fast explicit operator splitting 
method given above, and then condnct some nnmerical accuracy tests to verify 
our results. 



3.1 Error estimate 


We denote by u{x,y,t) the splitting solution satisfying exactly the scheme (2.3), 
and write t/j). := u{xj,yk,tn), '■= u{-, -^tn) and := u{xj,yk,tn)- We denote 
by and the discrete approximations of the operators Sj\f and Sc, respec¬ 
tively, and by the numerical approximation of satisfying 

Dehning a sample operator —)■ as I^u = {u{xj,yk))jk, we have 

Un _ where = {u E L^{Q) \ u is hl-periodic}. For the simple nota¬ 

tions, we omit the | or r following the symbols Sj\f, Sc, S^, or S^ below. 

To estimate the error, we need some lemmas. For the simplicity, we write Sj^u 
to mean and Sj^j-v to mean though the operators Sj,/ and S^ are 

actually nonlinear. We restate the accuracy of obtained in Section 2.2. 

Lemma 3.1. Under the condition (2.8), there exists a positive constant Ci, inde¬ 
pendent on T and h, such that 

WI^Smu - S^^l’^uW < C'ir(r=* + h^), Vm G 

We can derive the stability of S^ by using the result of Sj\f. 

Lemma 3.2. Given m E N. Under the condition (2.8), there exists a positive 
constant C 2 , independent on r and h, such that 

||5> - S^w\\ < ||n - wll + 2Cit{t^ + h^) + C 2 h^, Vn, w E 

Proof. Let v be some function, belonging to FT™,. (12), such that I^v = v, for 
example, the two-dimensional trigonometric interpolation of v in 12. Similarly, let 
w E Fr™j,(12) such that I^w = w. Using Lemma 3.1, we obtain 

|| 5 > - 5 >|| < || 5 > - I^S^fvW + Wl'^Sj^v - I'^Sj^wW + Wl'^Sj^w - 5 >|| 

< -U h‘^) + \\I^{Sj>/v - 5ArF;)|| -1 + h‘^). 

Since the L^-norm of an 12-periodic function on 12 can be approximated by the 
discrete L^-norm with spectral accuracy [26], using Proposition 2.1, we have 

||/"(5ArT-5ArFJ)|| < \\SMV-SMw\\L^(n) + Ch^ 

< ||F - w\\mn) + Ch^ < lln - m;|| + Uah™. 

Therefore, we obtain 

\\Sjfv - < ||n - m;|| + 2Cit{t^ + h^) + 

which completes the proof. □ 

Remark. In [14], the authors have proved the regularity of the solutions to the 
MBE equation (1.1) using the standard technique of Galerkin approximations. It 
says that u{t) E iF™j,(^) ^ > 0 if u{t)) E Fr™,.(12). With the similar proof, 

we can obtain Sj\f{t)u E iF™,,(12) and Sc(t)u E for any t > 0 provided 

u E iF™,.(12). Here we omit the detailed proofs and just use the results directly 
above. 
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The stability inequality (2.15) can be rewritten in the following form. 
Lemma 3.3. < e^||n||, \/v E 

The error estimate of the operator dehned in Section 2.3 can be proved in 
the framework of spectral method. 

Lemma 3.4. Given m E N and m > 1. There exists a positive constant C^, 
independent on r and h, such that 

Wl^Scu - < C'sIwUe^h™, Vm e 

Proof. We use the notations w{x,y,t), w{x,y,t) and w{x,y,t) representing 

OO OO OO N 

w{x,y,t)= w{x,y,t)= Upq{t)e-x‘^P^+^y\ 

p=—oo q=—oo p=—oo q=—N 

N N N N 

w{x,y,t)= w{x,y,t) = 

p=—Nq=—N p=—Nq=—N 

where Upq(t) and Upq(t) are given by (2.10) and (2.12), respectively, Assuming that 
w{-, •, 0) = w{-, •, 0) = w{-, •, 0) = w{-, •, 0) = M, we know that 


2N 2N 

Wl^Scu - S^I^uf = \w{xj,yk,T) - w{xj,yk,T)\‘^ < 4Ai + 4^2 + 2B, 

1=1 k=l 


where 

2N 2N 

Ai = h^^^\w{xj,yk,T) -w{xj,yk,T)\‘^, 
1=1 k=i 
2N 2N 

A 2 = h^'^'^jw(xj,yk,T) -w(xj,yk,T)j^, 

1=1 k=l 
2N 2N 

B = h'^^^\w{xj,yk,T) -w{xj,yk,T)\‘^. 
1=1 k=l 


We hrst estimate the terms Ai and A 2 . Since 

2N 2N , OO 


^. = ^^EE E E 

j=l k=l p=—oo |g|>W 

OO 

= 

p=—oo |ij|>Af 


Upqir)ex(P^^+<}yA 


< AL'^e^s ^ 

p=—oc |(j|>A'' 


U. 


pq 


p=—oo 


< AL'^e^s I iV-2”^ Y 

\q\>N 


u 


pq 
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and, similarly, 


= 4L^e6/V-=-” ^ ^ 9"”|S„(0)r, 

p=-oo |q|>Ar 


2N 2N N 

"*2 ='“'EE E E 

i=l k=l q=-N \p\yN 




N 


< 4 L^e 2 «iV E E 




q=-N |p|>Ar 


we obtain 


Ai + A2 < AL^e^^N E E -- 


N 


^pq 


p=—oo |g|>Af 
00 00 


+ E E p""i“" 

q=-N |p|>Ar 


<4L2e^iV-2™ ^ ^ 9nM2 




p=—oo q=—oo 
00 00 


< AL^ei-sN-^^ Y1 J2^P^ + (!^r\% 


^pq 


p=—oo q=—oo 


= AL^e^^N-^^\u\l^, 

where | ■ \m represents the semi-norm of 

We next estimate the term B. It is easy to obtain 


2N 2N N N 

^ S S (“P9(^) “ Upq{T))e'i^^^+'iA 

j=l k=l p=—Nq=—N 
N N 

= AL^ X] X] ~ 

p=—N q=—N 
N N 


\Upq[T) Upq[T) 


<AL‘^e2s 

p=—N q=—N 


U 


pq 


Now we look for the upper bound of B via the following fourth steps, 
(i) Magnify the sum 


N N 

^ E E 

p=—N q=—N 


^pq ^pq I 


A direct calculation leads to 

N / N-1 


^’“E E 

p=—N ' q=—N+l 
N , N-1 

s E ( E 

p=—N ^q=—N+l 


Upq Upq I , ^ 


+ ^ 'y ^ \2Upq 2u. 

q=±N 


pql 


~ ,2 1 

Upq Upq\ + 2 


\^Pq ‘^Upql + 2 \Upq\ 

q=±N ^ q=±N 
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< 


< 





\u. 


'pq\ 


q=±N 
. N 


nE E 


U 


pq\ 


p=-N q=±N 


E 


u 


pq\ 


N N 

= E E 

p=—N q=—N 

=: Di + D2 + -D3, 


p=±N 

N 


N 


+ E E 


u 


pql 


p=-N q=±N 
. N 


l^pq CpCqtlpq \ 


nE Em^+^E E 


\u. 


pq\ 


q=-N p=±N 


p=-N q=±N 


where Cp and Cq are defined as (2.14). 

(ii) Estimate the term D 2 + D^. Since 


1 


N 


N 


2 

1 


u 


pq\ ^ 


^’-^E E 

q=-N p=±N 
N 

E E p 

q=-N |p|>Ar 


q=-N 


^9 E 


2m\C: 12 


U 


-pql 


H>iV 

00 00 


2m|:rr |2 


u 


-pql ^ 


< 5^-"” E E p 


2m|;rs 12 


u 


-pql ) 


and, similarly, 


we obtain 


00 00 


D, < 5 ; 5 ; q 


q=—OQ p=—oo 


2m\Cr 12 


u 


-pql ) 


p =—00 q=—oo 


^ 00 00 ^ 

D2 + D3 < E E (P' + 

q=—oo p=—oo 

(iii) To estimate the term Di, we first prove that 

CpCqUpq ^pq T ^ ^ '^p-\-2rN,q-\-2sN • 

r^+s‘^f^0 

In fact, snbstitnting (2.9) into (2.12), we have 

2N 2N / CO CO 


CpCqUpq ^^2 


EE E E 


j=l k=l ^ r=—oo s=—oo 
2N 2N 00 00 

j=l k=l r=—co s=—co 

00 00 2N 2N 

r=—co s=—co j=l k=l 

CO CO 


4iV2 

1 


E E '^p+2rN,q+2sN 


r=—co s=—co 


(3.1) 


12 



= u 


■pq 


+ E 'Up+2rN,q+2sN j 


r^+s^^O 


Since 


2N 


r-p = 2lN, 




i=i 


0, r — p ^ 21N, 


here / is an integer. 


(iv) Estimate the term Di. Using the formula (3.1) and the Cauchy-Schwarz 
inequality, we have 

N N 

Di = E E E Up-\-2r N ,q+2sN 

p=—Nq=—N 

< E E |( E [{P + 2rNf+ {q + 2sNf\ 

p=—Nq=—N ^ h j. 2_|_^2 ^q 

^ ^ [(p + 2riV)^ + (g + 2siV)^] \Up-^- 2 rN,q+ 2 sN\‘ 


< max d E [(p + 2riV)^ + (g + 2siV)^] ”^Y 
bl,kl<^ Vo ,_o2^n / 


N N 

E! E/ E/ [(p + 2riV)^ + (g + 2siV)^] \Up+2rN,q+2sN\' 

p=—N q=—N r^+s'^^0 


< r E [(2^^ “ 

Vj,2^^2^Q / 

= 2N-^"‘Ml E 


■2|m 


|2 

Im 


r2+52^Q 


[(2r - 1)2 + (2s - 1)2]™' 


The series 

E 


< 


E 


1 


..+,V0 [(2-' - 1)" + (2" - l)d’‘ " ...fl.. (2’' - 1)"” + (2s - 1)"” 


1 m — 

- 2 ^ \2r - lM2s - 11 

r^+s^^O ' ' ' ' 


< oo, if m > 1, 


so we obtain Di < 2SN where S is the sum of the series above. 

As a result of (i)-(iv), we obtain 


B < 2{1 + AS)Lh^s N 


Ar-2mu|2 


So we obtain 


Wl^Scu - < 4(5 + 4A)L>|^eSiV-2"^, 

which leads to the expected result. □ 

Now we write the discrete L2-error estimate as the following theorem. 
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Theorem 3.1. Assume that Uq G with m > 1 and the condition (2.8) 

holds. If we set = I^uo, then the discrete If-error at T = nr is 

\\U^-u^\\<c(t‘^+ h^ +—y (3.2) 

Furthermore, if m > 6 and r ~ h'^, then 

\\U^-u^\\<C{t^+ h^). (3.3) 

Proof. Assume that u^~^, the numerical solution at t„_i-level, is given, then the 
discrete L^-error at t„-level should be 

11^7” - u^W < \\u^ - + \\U^ - u^\\. (3.4) 


The Strang splitting scheme (2.3) is second-order [24], which means that 

\\U^-U^\\ < CqT^ 


The second term in the RHS of (3.4) can be bounded as follows: 

I|t7- - u^W = Wl’^ScS^Scu^-^ - 

< WI^ScSmScu'^-^ - Sfl^S^Scu^-^W + WS^I^SmScu^-^ - 

< C^\SMScu^-^\mei-^h^ + ei-^Wl'^SuScu'^-^ - (3.5) 


where the last inequality is the consequences of Lemmas 3.3 and 3.4. Besides, 




< WI’^SmScu^-^ - Sjfl’^Scu^-^ + WSl^l’^Scu^-^ - 

= + WSj^il'^Scu^-^) - 

< 3 Cit{t^ + h^) + Cah™ + Wl’^Scu'"-^ - S’fu^-^\\, 

(3.6) 


where the last ineqnality is the consequences of Lemmas 3.1 and 3.2. Fnrthermore, 

Wl'^Scu"-^ - S'fu^-^W < \\I^Scu"~" - S’fl^u^-^W + - S'fu^-^W 

= \\{I^Sc - S^I^)F^-^\\ + \\Sf{U^-^ - n"-')|| 

< C'3|n"“Ve®/i" + e^||R"“^ (3.7) 


where we use the fact U"' ^ Combining (3.5)-(3.7) with (3.4), we obtain 

\\U" - u"\\ < e* + 3C'ie*r(r3 + h^) + (Ca + C'3 C't( 1 + e*))e^h’", 

where Ct = m.eex{\Sj^Scu^\m-, \u^\m : 0 < /c < n}. Setting = \\U'^ — u^\\ and 

hm 

G = SCie^ir^ + h^) + {C 2 + CsCril + e*))e*—, 

we have 

F^ <e^F^-^ + tG, n = l,2,.... 
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Using the Gronwall’s lemma and the fact e^ — 1 > a; (a; > 0), we obtain 

T 

pn < ^^po ^ - )_Q < ^ 25e^sG. 

e^s — I 

Since F° = \\U^ — m°|| = 0, we obtain 

||U” - All < 2Se^ (^3C'ie^(A + A) + (A + C'3 C't( 1 + e^))e®—j, 
and thus 

. um. 

WU'^ - All < CoT^ + 2(5eS (^3C'ie®(A + A) + (A + C'3C'r(l + e®))e®—j. 

If r < min{4(51n2,1}, then e® < 2, so we obtain 

T T T 

\\U^ - All < (Go + 12(5e^Gi)A + 12(5e^GiA + 4(5e^(G2 + SCsCt) — , 

T 

which implies the estimate (3.2). 

Furthermore, we set the step r ~ to obtain 

UTn 

^ Lm—2 1 

— ~ /r ~ 7-2 

T 

As long as m > 6 holds, we obtain the error estimate (3.3). □ 

3.2 Accuracy tests 

Now we carry out the accuracy tests on the equation (1.1) with 5 = 0.1, T = 1, 
n = (0,27r) X (0,27r), and 

Uf){x,y) = 0.1 (sin 3a; sin 2y + sin 5a; sin by ), 

which is a classical example studied either theoretically or numerically [14, 17, 19, 
29]. We take the numerical solution obtained with r = 5 x 10“® and J = 2048 
as the “exact” solution. The tests are conducted with different spatial scales, and 
the time step is set to be r = where the constant Cq is chosen to render 

r = 0.005 when J = 128. Table 3.1 shows the discrete L^-errors implying the 
accuracy nearly 0{t‘^ + A), which is consistent with Theorem 3.1. 


Table 3.1; The discrete L^-errors with different spatial scales and time steps. 


J 

r 

l|e(AII 

e(J/2)| 

l|e(J)|| 

log2 - 

\\e{J)\\ 

128 

5 X 10-3 

1.0278 X 10-5 

* 

* 

256 

1.25 X 10-3 

9.5361 X 10-^ 

10.7779 

3.4300 

512 

3.125 X 10-^ 

6.5869 X 10-3 

14.4774 

3.8557 

1024 

7.8125 X 10-5 

2.4026 X 10-*^ 

27.4156 

4.7769 
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4 Numerical experiments 


Example 4.1. We consider the one-dimensional MBE model 


^'^XXXX’! 

(x,f) e (0,12) X (0,T], 

is 12-periodic, 

t G [0,T], 

f irx 2'xx \ 

X e [0,12]. 

M(a;, 0) = 0.1 ( sm - h sm-h sm vrx 1, 

\ 2 3 / 


The evolution of this initial-boundary-value problem is studied theoretically via 
the perturbation analysis [14] to observe the morphological instability due to the 
nonlinear interaction. It is also a classical example for the numerical experiments 
in the case 5 = 1. Here we will present the results obtained by the operator 
splitting method given in Section 2. Fig. 4.1 shows the results of the case 5 = 1 
with J = 128 and r = 0.1, which is consistent with the existing work [14]. 



(a) height: t = 0 

0.4 I-^^^^^- 

0.2 

- 0.2 

-0.4 I-^^^^^- 

0 2 4 6 8 10 12 

X 

(d) height: t = 20 



(g) height: t = 100 


0.41- ^^^^^ -1 

0.2 

0 - 

- 0.2 

-0.4 1 - ^^^^^- 

0 2 4 6 8 10 12 


(b) height: t = 0.5 



(e) height: t = 30 



(h) gradient: t = 100 


0.41-^^^-1 

0.2 

0 - 

- 0.2 

-0.4 1 - ^^^^- 

0 2 4 6 8 10 12 

X 

(c) height: t = 15 



(f) height: t = 60 



(i) energy: 0 < t < 100 


Figure 4.1: Example 4.1: The results of the case 5 = 1. 


Besides, we also present some results from reducing 5 to 0.1, 0.01, and 0.001, 
respectively. The results are summarized in Figs. 4.2-4.4. The hrst plot in each 
hgure presents the height u{x,t) at some time t, the second one shows the cor¬ 
responding gradient Ux{x,t), and the third one plots the evolution of the energy 
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1.5 



X X t 

(a) height: t = 200 (b) gradient: t = 200 (c) energy: 0 < t < 200 

Figure 4.2: Example 4.1: The results of the case 6 = 0.1 obtained with (J, r) = 
(128,0.01) (solid line) and (J, r) = (256,0.005) (dash line). 



(a) height: t = 500 



X 


(b) gradient: t = 500 


3 

2.5 
2 

1.5 


0.51 -^^^^- 

0 100 200 300 400 500 

t 

(c) energy: 0 < t < 500 


Figure 4.3: Example 4.1: The results of the case 6 = 0.01 obtained with (J, r) = 
(256,0.001) (solid line) and (J, r) = (512,0.0005) (dash line). 



(a) height: t = 1000 (b) gradient: t = 1000 (c) energy: 0 < t < 1000 


Figure 4.4: Example 4.1: The results of the case 6 = 0.001 obtained with (J, r) = 
(512,0.0001) (solid line) and (J, r) = (1024,0.00005) (dash line). 


Fig. 4.2 presents the results of the case 6 = 0.1 with (J, r) = (128,0.01) and 
(J, r) = (256,0.005). We find that both solutions have few differences between 
each other, so we are convinced that the results presented here are credible. The 
energy decreases hardly after t = 200 so that we view the solution at t = 200 as 
the steady state. It is observed from the left and middle graphs that there are 
two complete waves in the steady state whose gradients do not exceed the range 
between —1 and 1. 

Fig. 4.3 gives the results of the case 6 = 0.01 with (J, r) = (256, 0.001) and 
(J, t) = (512, 0.0005). Fig. 4.4 gives the results of the case 6 = 0.001 with (J, r) = 
(512,0.0001) and (J, r) = (1024,0.00005). Likewise, we can trust these results. 
The solutions at the steady states present more waves in the considered domain 
than those above, and the gradients still locate in the interval [—1,1]. 

From the gradient graphs of Figs. 4.1-4.4, we find that the smaller 6 is, the 
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more points on the gradient curves locate at the horizon lines y = 1 oi y = —1. 
This is a consequence of the competitions between the Ehrlich-Schwoebel effect 
and the dissipation mechanism of the energy functional 



It is seen that the Ehrlich-Schwoebel effect selects the slope \ux\ = 1 while the dis¬ 
sipation term weakens the selection. To reduce 6 means to weaken the dissipation 
effect, or equivalently, to strengthen the slope selection. And thus, the gradient 
interfaces connecting —1 to 1 turn steep and the solution curves turn sharp there. 

Example 4.2 (Coarsening dynamics). We simulate the two-dimensional MBE 
model (1.1) with 5 = 0.1 on the domain Vt = (0,100) x (0,100) with 512 x 512 
grid. We set the initial data as a stochastic state by a random number varying 
from —0.001 to 0.001 on each grid point. The time step is set to be t = 0.01. 


This example is aimed to verify the power laws for the energy evolution and 
the height growth. 

Fig. 4.5 presents the contour lines of the free energy 


Ffree:=^(|Vwp-l)2 + ^|Aw 

at t = 100, 2000, and 30000. 



Figure 4.5: Example 4.2: Contour lines of the free energy Ftee- 


Fig. 4.6 shows the evolution of the energy and the interface height. The (a) 
presents the power law for the evolution of the energy. The energy curve is plotted 
in log-to-log scale and nearly parallels to the dash line fa, which suggests that 
the energy evolves in time as the power law Cf^ with a — 1. The (b) presents 
the power law for the growth of the interface height h{t), which is dehned by 

h{t) = j u‘^{x,y,t)dxdy 

Again, the height curve is plotted in log-to-log scale. The growth of the interface 
height approximately obeys the power law Ct^ with (3 ^ which is consistent 
with the existing works (see, e.g., [17, 19, 29]). 
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^-1/3 



- 1/3 law 


10 ^ 10 ® 10 "^ 


t 

(a) energy evolution 

Figure 4.6: Example 4.2: The power 
interface height. 



(b) height growth 


for the evolution of the energy and the 


5 Conclusions 

In this work, we investigate the error estimate of a fast explicit operator splitting 
method for a nonlinear fourth-order diffusion equation modeling epitaxial growth 
of thin films. The convergence order in discrete L^-norm is proved 

theoretically and verified numerically. For the nonlinear subproblem, we construct 
a 25-point center-difference scheme in space and use the third-order explicit SSP- 
RK scheme in time. Since fewer points are involved in the scheme at each node 
compared to the 33-point center-difference scheme presented in [1], the restriction 
for the stability may be reduced. We carry out several numerical experiments to 
verify the efficiency of the derived algorithm and present some results for small 
(5’s with the time step r = 5/10. Furthermore, we hnd numerically the coarsening 
exponents for the energy evolution and the height growth are —1/3 and 1/3, 
respectively. 
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